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Abstract 

A solution based grid adaptation method that combines elements of the multigrid method for solution acceleration 
and the domain decomposition philosophy for grid optimization is described. Unlike other solution based adaptive gridding 
schemes, wherein the overhead of recomputing the grid and re-evaluating the solution on the adapted grid leads to higher 
computational costs compared to a non-adapted calculation, the present methodology reduces the computational time 
required to obtain the solution. The computational effort involved in the present calculation is significantly lower than a 
non-adapted calculation that utilizes the multigrid method purely as a convergence acceleration tool. In addition to 
convergence acceleration, the multigrid framework provides a mechanism of information transfer from regions wherein grid 
refinement is specified to unrefined coarse grid regions. The basis for domain decomposition in the current procedure is the 
variation in grid refinement requirements for each coordinate direction in different portions of the flow field. The method 
is demonstrated herein on an efficient set of governing equations termed the reduced Navier Stokes equations, applied in 
conjunction with a set of physical boundary conditions. The governing equations are discretized through a pressure based 
flux splitting procedure that is uniformly applicable from incompressible to supersonic Mach numbers. 

1. Introduction 

Although significant progress has been made in the field of high speed computing there is still a need for more 
efficient algorithms to solve large scale fluid dynamics problems that require considerable amounts of computer storage and 
computational time. Complex three dimensional flow computations and direct numerical simulation of large Reynolds 
number flows are examples of such problems. Strong viscous-inviscid interactions, associated with turbulent or high 
Reynolds number (Re) laminar flows, are quite frequently characterized by the appearance of large flow gradients that are 
most significant in small or 'thin* domains of finite extent, and in one or more directions, e.g., boundary or shear 
layers/regions, triple deck structures and vortical or recirculation zones. Outside of these regions, the flow is generally 
diffused or inviscid and the flow gradients are less severe. However, the flow character in these 'smoother' regions, which 
generally encompass the major portion of the flow domain, can be significantly influenced by interaction associated with the 
thin high gradient viscous layers. In this paper we propose an approach that optimizes the overall calculation from several 
standpoints and is applicable to a large class of flow problems ranging from incompressible to supersonic flow regimes. An 
efficient set of governing equations that effectively represent the asymptotic behavior of the full Navier-Stokes equations are 
applied. The equations referred to henceforth as the reduced Navier-Stokes (RNS) equations are discretized through a finite 
difference scheme that is derived through a pressure based flux splitting approach. A general solution based grid adaptation 
methodology that also ensures grid convergence efficiently is described. This approach combines the efficiency of the 
multigrid method and a domain decomposition technique that provides a mechanism to split the overall flow domain into 
subdomains within which adaptive grid prescription can be specified based on the local behavior of the flow field. The 
computational time required for most adaptive grid schemes is typically larger than that required for a non-adaptive 
calculation. This is not the case in the present approach. The adapted grid computation requires less time than a non- 
adapted calculation that applies the multigrid method for convergence acceleration. The multigrid method provides 
effective communication between the disparate flow domains through global coarser grids, and at the same time maintains 
all conservation properties. 

Many approximations of the full Navier-Stokes equations have been developed over the years, see for example, 
Davis et. al. (Ref. 1). In the present study the Navier-Stokes equations are represented by an implicit lowest-order reduced 
Navier-Stokes (RNS) system and a purely diffusive, higher-order, deferred-corrector. The equations are discretized through 
a pressure based flux split finite difference approach. A trapezoidal or box-like form of discretization ensures that all mass 
conservation properties are satisfied at interfacial and outflow boundaries, even for a primitive-variable, non-staggered grid 
computation. In the past, the RNS equations have been used to solve a variety of flow problems, for example Rubin et. al. 
(Ref. 2, 3), Rosenbaum et. al. (Ref. 4), etc., have applied these equations in conjunction with a global pressure relaxation 
solver for a range of low speed incompressible to low Mach number subsonic flows, for both internal and external 
geometries. Ramakrishnan et. al. (Ref. 5), Lai et. al. (Ref. 6), etc., have extended this methodology to unsteady 



compressible flows. Pordal et. al. (Ref. 7) have computed time accurate flow solutions, including supersonic inlet unstart, 
using the RNS system and methodology. 

The discretized equations are solved through a global pressure relaxation approach as described in Rubin et. al. 
(Ref. 3). This is a space marching approach wherein at each (approximately) streamwise location the discrete, non-linear 
equations are solved exactly based on known values (or current iterate) at the previous spatial station and guessed (or 
previous iterate) values at the next streamwise location. As with most iterative methods the convergence rate of this 
multisweep procedure, that improves the pressure field in each iteration, slows down significantly as the number of mesh 
points especially in the streamwise or marching direction is increased. Himansu and Rubin (Ref. 8, 9) have applied the 
multigrid method in a semi-coarsening mode to significantly accelerate convergence for a wide range of two and three 
dimensional external flow problems. In this approach the grid is coarsened only in one direction when proceeding from one 
multigrid level to the other. This approach was found to be more effective especially for turbulent flows and laminar flows 
with strong viscous-inviscid interaction and flow separation. The inability of the coarse grids to capture these flow features 
make the coarse grid corrections to the fine grid solution ineffective and the multigrid procedure breaks down. Some of 
these effects can be attributed to the non-uniformity of the grids typically used in the direction normal to the body surface. 
In fact, two dimensional laminar flow examples with strong separation and viscous inviscid interaction have been 
successfully computed through the present adaptive gridding multigrid domain decomposition approach wherein coarsening 
in both directions have been applied successfully. In this approach only uniform grids were used in all grid levels. In the 
current investigation, this grid adaptation procedure has been extended to three dimensional flows. For laminar flow 
computations grid coarsening/refinement is allowed in all three directions. For turbulent flow computations the semi- 
coarsening muitigrid approach (with the adaptive element) was found to be most effective because use of uniform grids can 
be uneconomical in three dimensions even though adaptive refinement is prescribed. Significant gains in computational 
costs are achieved, over single grid calculations, with either a full coarsening or semi-coarsening multigrid technique. 

Grid adaptation approaches can be broadly classified into two basic approaches: (1) Dynamic redistribution of a 
preset number of grid points by specifying a solution based criterion that minimizes the discretization error, and (2) the 
addition of grid points in regions where the truncation error is large. Each of these techniques has advantages and 
disadvantages. The first approach guarantees better accuracy; although, the overall computational time may not necessarily 
be lower than that for a non-adapted calculation. In addition, the level of accuracy is constrained by the number of grid 
points used in the initial grid. The additional computational cost associated with computing the grid distribution is high 
(Thompson et. al., Ref. 10). Also, special care must be taken to ensure a compromise between conflicting requirements of 
orthogonality, smoothness or lack of skewness, and concentration. A variational approach that couples the solution of the 
flow variables with that of the grid is by far the best methodology when dealing with moving grids. The second approach 
of overlaying grids also has several variants. In the CHIMERA method, first developed by (Steger et.al.. Ref. 1 1) a major 
grid is generated about a main component of the configuration and minor grids that may be generated through different 
transformations, are overlaid. This method provides the flexibility of handling complex, multi-component geometries by 
using simple grid generation tools. It also allows for the possibility of local refinement in regions of large gradients. The 
absence of a global coarse grid can hinder convergence of the solution, in particular, due to inconsistencies in grid accuracy 
of the different patches. Other overlaid mesh systems have been used previously by a number of researchers. Magnus et. al. 
(Ref. 12), used overlaid grids to achieve improved computational accuracy for solving transonic flow about airfoils. Berger 
et al. (Ref. 13) has developed a numerical procedure which automatically inserts overlaid grids to resolve high gradient 
regions to two dimensional convection problems. With overlaid grids, conservation properties at grid interfaces can be of 
major concern. 

Local refinement procedures necessitate proper treatment of grid interfaces to ensure mass conservation. 
Typically, boundary conditions on the local grids are obtained by interpolating the solution obtained on coarser grids. In 
full Navier-Stokes calculations boundary conditions are required on all three velocity components at non-wall boundaries. 
Hence simple linear interpolation schemes are not sufficient to ensure mass conservation within local sub domains. In the 
current RNS methodology, only two of the three components of velocity need to be prescribed at a free boundary. The third 
component is computed from the RNS system with a trapezoidal discretization of the continuity equation. 

Using the multigrid method as a framework for adaptive grid refinement was first suggested by Brandt (Ref. 14). 
The development of this methodology is very logical, since the multigrid method already comprises a hierarchy of grids 
ranging from coarse to fine. Conventionally, the multigrid method is applied to accelerate the convergence of a discretized 
elliptic problem on a given mesh. This framework can be conveniently adapted to provide optimal refinement, by limiting 
the extent of the finer grids to regions where the error, or lack of grid convergence thereof, is still large. Brandt refers to 
this as ‘segmental’ refinement. This ideology forms the basis of the current work. Earlier research, on similar approaches, 
has been presented by Fuchs (Ref. 15) and Thompson et. al. (Ref. 16): although, in their approaches no attempt is made to 
distinguish between refinement requirements in the different coordinate directions. This has been achieved by the authors in 
previous papers (Refs. 17,18) for a variety of two and three dimensional flows, that are primarily represented through a 
Cartesian coordinate system. The need for grid refinement in each coordinate direction is tracked through truncation error 
estimates of specific derivatives. 
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Adaptive gridding requires some form of error estimation. A number of different approaches have been applied 
by previous researchers, for the purpose of grid adaptation, as well as, for a means of estimating the quality of the results. 
Typically most error estimation or grid convergence checks require a solution on two different grids with different 
resolution. Brandt (Ref. 14) recommends estimate of the coarse grid truncation error as a refinement criterion. This can 
easily be estimated within the multigrid algorithm. However, this estimate does not provide information about which 
coordinate direction lacks refinement. In order to obtain such information, the test would have to be performed multiple 
times; each time the finer of the two grids has to be refined in just one coordinate direction. This again adds to the overhead 
of the calculation. It also makes it impractical as a criterion for adaptive refinement, where the primary goal, apart from 
better accuracy, is to reduce the overall computational cost. Hence in the present study, a simple and inexpensive procedure 
based on Richardson extrapolation is applied to estimate truncation errors of key derivatives. These are then used as criteria 
in the adaptive refinement process. 

Domain decomposition can be viewed as a divide and conquer philosophy, where sub problems corresponding to 
each sub domain can be solved using any appropriate efficient solver. There are interactions between sub domains, that 
involve passage of suitable information through an iterative procedure. It is recognized that when there are several sub 
domains, the presence of a global coarser grid can greatly enhance the rate of convergence of the iterative procedure. Since 
this is the underlying philosophy of multigrid methods, the current segmented multigrid domain decomposition concept 
weaves these two ideas to achieve significant gains in computational time over conventional multigrid procedures. The 
basis for domain decomposition in the present study is disparity in refinement requirements in different portions of the flow 
field. 

In the present study, results are presented for several three dimensional internal flow configurations for both low 
and high speed regimes. The discretization scheme is uniformly valid for the entire Mach number range from 
incompressible to supersonic flows. 

2. Solution Procedure 

In this section a detailed description of the adaptive gridding procedure is presented. Details relating to the grid 
structure, truncation error estimation and multigrid implementation are described. In addition some aspects of the data 
structure necessary for enhanced implementation of the local refinement procedure are described. 

2.1 Grid Structure 

The overall solution algorithm comprises several grid levels that form part of the multigrid hierarchy. The calculation 
starts on a fairly coarse grid for most problems. In calculations wherein the multigrid method is applied purely as a 
convergence acceleration tool, the finer multigrid levels comprises a single global grid that is determined apriori. The aim 
of the present approach is to optimize the computational grid. Hence, extent of the finest grid is not known apriori 
necessitating the calculation to start on coarsest grid levels and proceeds to finer levels in the multigrid hierarchy. The finer 
multigrid levels are generated through the adaptive refinement process. Assume grid adaptation is initiated at level k 
(typically k = 2) in the multigrid hierarchy; level I being the coarsest. The grid starting at level k can comprise series of 
subdomains each refined with respect to the k-1 grid level in one or more coordinate directions. The subdomains within a 
multigrid domain could be abutting each other or disjoint. This flexibility allows disparity in refinement requirements in 
different regions of the flow field to be handled efficiently. 

The choice of the coarsest grid has to be specified with care and is determined entirely by the complexity of the 
geometry and the associated flow field. In fact, some calculations have been performed with the first mesh point, from wall 
boundaries, lying outside of the thin viscous or boundary layer (Ref. 17-19). Unlike laminar flows, in complex three 
dimensional turbulent flow calculations it may not always be possible to obtain solutions on very coarse grids, especially in 
directions normal to the walls when wall functions are not applied which is the case in all the computations presented herein. 
Similarly for flows involving shock waves, excessively coarse grids cannot be used in the streamwise direction. The 
coarsest grid should not contain a large number of grid points, as this will result in high overhead cost for the truncation 
error calculation, and ajess than optimal application of the adaptive gridding approach. These considerations are primary in 
the choice of the coarsest grid distribution. 

The first two grid levels necessarily span the entire computational domain. A relatively coarse mesh size is used in all 
coordinate directions, along which adaptive refinement is to be performed. Thus if adaptive refinement is to be performed 
only along the ^-direction, then a coarse grid is prescribed in that direction. In the other two directions, there is the option 
of using coarse or fine grids. This depends on the choice of the multigrid algorithm to be applied, i.e., (i) full coarsening 
multigrid or (ii) semi-coarsening multigrid. In problems where coarse grids are unallowable due to the flow physics 
considerations mentioned earlier, the semi-coarsening multigrid philosophy appears to be a better choice, although full 
coarsening is usually possible for one or two levels. Four or five full coarsening levels are possible in laminar flows and 
turbulent flows in simple geometries. 

The global pressure relaxation procedure is applied to obtain the solution on the first two coarse grids. This 
constitutes less than 5-10 % of the overall computational cost. Based on the truncation error estimation process to be 


3 



described in the next section, the next multigrid level is determined. In general, level k in the multigrid hierarchy is 
determined based on the truncation error estimation using solution on level k-I and k-2. After the extent and grid structure 
of level k is determined, the calculation proceeds sequentially from one subdomain to another. One relaxation sweep on a 
fine, adapted multigrid level comprises of one sweep in each of the subdomains in that level. The decision to transfer the 
calculation to another grid level depends on several factors. These will be discussed in detail in the section on multigrid 
implementation. Upon convergence of the current finest multigrid level, the truncation estimation process is repeated with 
the solution on the two finest grid levels and the domain structure of the next multigrid level is defined. Since, truncation 
error is only determined for all points in the finest grid, the regions refined will naturally be a subregion of the finest grid 
although there could be further segmentation based on the local truncation error in each coordinate direction. It is also 
possible that the entire region occupied by level k is refined in level k+1. 

2.2 Truncation error estimation and refinement strategy 

Error estimation is an essential element of all adaptive gridding methods. Zeeuw et. al. (Ref. 20) adopt a criterion 
based on the change in the solution and they fix the percentage of grid points that should be refined. Thompson et. al. (Ref. 
16) apply a criterion based on the global truncation error of the system of equations. Brandt (Ref. 14) has shown that given 
the solution on two grids, the global truncation error of the coarse grid can be approximated with the fine grid accuracy. 
This type of truncation error evaluation provides a reasonable guideline to ensure grid convergence, but it does not provide 
information about which coordinate contributes significantly to the truncation error. 

In the present study one of the primary aims is to optimize the grid, in each coordinate direction independently, 
wherever possible. Information relating to the degree to which a given coordinate gradient contributes to the overall 
truncation error is required. Since this is not easily achieved through a global truncation error procedure, error estimates in 
certain key derivatives are evaluated and applied as the criteria of refinement for each coordinate direction. For the 
streamwise (E) direction, the truncation error in the p^ term is used as the criterion for both compressible and 
incompressible flows. Most of the flow phenomena to be considered herein, e.g., separation, reattachment, shocks, strong 
geometric curvature, etc., are associated with large pressure gradients. Thus this term can be an effective mechanism in 
identifying such regions. For the normal (rj) direction, the truncation error in vorticity or u^ is used as the criterion for 
local refinement. Boundary layers are usually associated with strong vorticity gradients and u^ is the most significant 
contributor. Additional gradient parameters can be added when necessary, e.g., u^ , if ^ adaptation is desired. New 


criteria can be defined to satisfy any specific application. For example, if it is desired to apply the adaptive refinement 
philosophy for predicting laminar instability and transition, a criterion based on the frequency of disturbance or the rate of 
growth of the instability can be applied. 

The truncation error in the ‘key’ derivatives is obtained with a Richardson extrapolation procedure. This provides a 
simple and reliable estimate of the truncation error for specific derivatives. Conventionally, Richardson extrapolation is 
applied to improve order of accuracy. An appropriate combination of solutions is obtained on two different grid levels and 
the leading truncation error terms are eliminated. In the current approach, the same philosophy is applied to calculate the 
leading truncation error term. For example, consider the p^ term. The discretization of this term, is given by (Ref. 17, 17). 
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The simplifying assumption o)j _ i/2 = o ) i+ | /2 has been made in equation (1). HOT refers to higher order terms. A similar 
expression for the coarse grid, with grid spacing 2A§, is 
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The superscript ‘c’ and ‘f in equations 1 and 4 refer to the fine and coarse grids respectively and the subscript I refers to 
the same physical grid location as i on the fine grid. Letting T p ^ = and subtracting equation (1) from (4) we 


obtain. 
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A similar expression can be derived for the u^ truncation error. 

This approach provides a very reliable method of identifying regions where truncation error is large although the 
absolute value of truncation error calculated may not be accurate. In the present approach the aim is to restrict grid 
refinement to regions of large truncation error. Richardson extrapolation provides this information in a reliable fashion. 
The absolute values of the truncation error is not needed in the present adaptation procedure. 

The truncation error obtained through this procedure is normalized with the maximum value over all the grid points in 
the current multigrid level. The normalized truncation error always ranges between 0 and 1. This unifies and simplifies the 
setting of tolerances for the various truncation errors. The values of non-normal ized truncation errors can vary significantly 
from problem to problem and will be dependent on parameters such as the Reynolds number, grid spacing etc. The choice 
of the tolerance for the normalized truncation error is based on range of mean, maximum and minimum values of the 
truncation error field. If the truncation error field is fairly uniform then a lower tolerance must be chosen to ensure that not 
too many grid points are refined. This is usually not the case for complex viscous flows. Even for simple boundary layer 
flows the distribution of truncation error in the discretization of the gradients normal to the surface depict a large variation. 
Similarly, large truncation errors are associated with streamwise pressure gradient in regions of flow separation, 
reattachment and shocks waves. The variation in the truncation error in the pressure gradient is not as much as the 
truncation error in the normal viscous derivatives. Hence a higher tolerance (0. 1-0.3) is set for the normalized truncation 
error in the pressure gradient term whereas a much lower tolerance (0.01-0.1) is set for the normalized truncation error in 
the viscous terms. This ensures that too many grid points are not refined. Under-refining the grid will result in poor 
solutions and/or non-convergence of the finer grid levels. This is discussed in more detail discussed in Ref. 22. 

Coarse Fine 


Case 1 

Coarse Fine 



Case 2 


Case 3 



Figure 1 Sample refinement patterns. 

Figure 1 shows typical grid structures on two successive multigrid levels. In Case 1 , the grid is refined locally in the ^ 
direction. Similarly in Case 2. the grid is locally refined in the T| direction. The finer grid contains only one subdomain in 
both these cases. For Case 3 however. 5 subdomains are identified - two with T) refinement, one with ^ refinement, and 
two with refinement in both £, and r| directions. Note that in Case 3, there are both abutting and disjoint subdomains. 

The process of calculating the truncation error and defining new sub domains is repeated for subsequent multigrid 
levels. The treatment of each sub domain is the same as for the global domain. The type of refinement associated with a 
new sub domain is a subset of the type(s) of refinement associated with its parent subdomain. In each subsequent level the 
amount the sub-region where refinement is prescribed decreases. Each multigrid level is decomposed into a discrete 
number of abutting or disjoint sub domains, for which refinement is specified in one or more directions. This pattern is 
defined by the directional refinement criteria as determined by the truncation errors in the two key derivatives mentioned 
above. Each successive multigrid level derives part of its topology from the subdomaining pattern of the coarser 
predecessor. 

Although within individual subdomains the grid is structured, when all the multigrid levels are overlaid, the effective 
grid has the appearance of an unstructured grid. In the global pressure relaxation procedure a line/plane is solved implicitly 
in the marching direction. This requires each computational unit, in the present case each subdomain, to be a rectangular 
region in the (|,r|,Q computational system. During the truncation error estimation process, all grid points which need 
refinement in the | and T} directions are tagged based on individual refinement criterion. A marching sweep is performed 
in the £ direction and all (t|,Q planes that contain at least one point needing refinement is marked. These planes are then 
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grouped into individual blocks; some blocks require refinement and some do not. Once this blocking is completed, each 
block is searched for points needing rj refinement. During this process, blocks are created in the T| direction. A similar 
strategy is followed as the block creation process in the % direction. Within each block in the direction, intersections with 
T) blocks are determined and individual sub-domains that would form the next multigrid level are determined. In each 
subdomain appropriate combination of refinements are performed. During the block creation process in the £, and r| 
directions, a tolerance of two stations is provided before creating a new block i.e., a new block is created only if two 
successive stations are encountered without points needing refinement. No attempt is made to adaptively refine in the £ 
direction in the current work. This means that either full refinement is performed or fine grids are used in this direction 
even on coarse multigrid levels. During the process of grouping points needing refinement into rectangular regions, points 
which do not require refinement can get refined. As a result the final grid is not necessarily optimal but it is optimized. If 
a explicit solver is used to solve the system of equations, unlike the semi-implicit solver presented herein, the refinement 
process can be optimized further, although the efficiency of explicit solvers for complex flows with strong viscous inviscid 
interaction is debatable. 

Some guidelines are followed during the refinement process. Within each sub domain, refinement is specified based 
on the overall truncation error level for all sub domains of that multigrid level. The type of refinement performed within a 
new sub domain must be a subset of the type(s) of refinement performed within its parent. For example, further refinement 
(at a higher multigrid level) in the streamwise direction is not performed in a sub domain that was refined only in the 
normal direction in a coarser multigrid level. Refinement cannot be re-initiated in a region where it was terminated and the 
extent of any sub domain cannot exceed that of its parent. This strategy is not suited for unsteady flows where the region 
needing refinement has to vary as significant flow features are convected with time. A methodology of shifting the fine 
grid regions to capture the transient flow behavior has to be devised for such situations. A complete redefinition of the 
adapted grids can also be performed for each time step. 

Since the refinement is performed locally, ‘hanging nodes’ will appear to be present in the grid when the locally 
refined grid is overlaid on the global coarse mesh. As described in the boundary condition section to follow, the values at 
these interior boundaries are obtained by interpolation of either the coarse grid solution values or solution obtained on 
neighboring subdomains which may be sharing the interface in the same multigrid level. Although the nodes at the 
interface appear to be hanging nodes, within a single sub-domain there are no hanging nodes, since the grid is uniformly 
refined throughout the sub-domain. Thus, as the calculation proceeds on a given sub-domain, no special treatment is 
necessary. There is no difference in the implementation of the global pressure relaxation procedure for local sub-domains 
and global coarse grids. The refinement process does not impose any constraints to control aspect ratio of cells or ratio of 
refinement levels between adjacent grid points. Thus when the overall grid is viewed as a union of all the multigrid levels, 
there could be physically adjacent points that differ in refinement by more than one level. But when the calculation 
progresses, there are no sudden jumps in the grid because the level of refinement within each domain is uniform. The 
global pressure relaxation algorithm is applied sub-domain by sub-domain. The boundary conditions are updated only 
during visits to the coarser grid. Although, typically size of the sub-domain wherein refinement is prescribed decreases, 
this is not a required in the present method. Situations wherein the parent domain is refined to form the subsequent 
multigrid level is allowable. 

For geometries that require grid transformation, the finest grid is generated in all three directions. Thus the maximum 
number of allowable grid levels is prescribed apriori. Typically a grid that is extremely fine can be generated and stored. 
When a grid level is defined adaptively, with respect to its parent, then only a portion of this global fine grid corresponding 
to that level is required. By ensuring that all coarse grid points are subsets of the finest grid level, the geometric definition 
is preserved in coarser levels. Figure 2 shows the adapted grid obtained for a 90° bend duct. Note that the finest level 
consists of four separate sub domains. These comprise portions of a global grid, that is generated for the geometry. A 
simple algorithm to search for the starting locations of each sub domain and to read in the appropriate portion of the global 
grid has been implemented. This aspect of the methodology leads to a significant advantage especially for complex 
geometries wherein it is often difficult to limit the number of grid points and meet all the grid clustering features desirable 
to obtain a grid independent solution. The current adaptive gridding procedure takes away this burden by allowing a fine 
grid to be generated throughout the domain without the constraint of limiting the number of grid points. Portions of this 
global fine grid are then extracted to define the finer multigrid levels. The adaptive refinement process creates a grid that 
resembles an unstructured grid, but retains the advantages of a structured grid. 

2.3 Multigrid Implementation 

The multigrid technique is used here to solve a system of discrete equations that define an elliptic boundary value 
problem on a hierarchy of successively finer grids. This technique can be viewed in two different ways; (i) the coarser 
grids are correction grids that accelerate the convergence of the finest grid by efficiently removing low frequency error 
components, or (ii) the fine grids are correction grids that improve the accuracy of the coarser grids by introducing fine 
grid transfer functions to the coarse grid discrete system. Both view points are relevant to the present study. On coarse 
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grids levels, wherein local refinement has not yet been initiated, the component (i) is significant. Through component (ii), 
the finer grid levels improve the accuracy of the coarse grid solution, even outside of regions where refinement is 
performed. 



The multigrid technique has been applied to a wide variety of compressible flow problems by Jameson et. al. (Ref. 
23), among others. Acceleration of compressible and incompressible flow solution procedures have been presented by 
Shyy et. al. (Ref. 24). The technique has been applied to unstructured grids by Mavriplis (Ref. 25). Himansu et al., (Ref. 
8, 9) have applied the multigrid technique in a semi-coarsening mode to accelerate the convergence of the global pressure 
relaxation procedure for a wide range of incompressible and compressible external flows in both two and three dimensions. 
In this approach, the grid is coarsened only in the £, direction when transferring from fine to coarse grids. The details of the 
implementation of this multigrid technique for the global pressure relaxation procedure can be found in Himansu et. al. 
(Ref. 8, 9) and in Srinivasan (Ref. 19). Certain elements are reviewed here for the sake of completeness, and to illustrate 
some interesting findings with regard to the multigrid restriction operation. The discrete equation that is solved on portions 

a k-l 


of the coarse grid, that have been refined is of the form 


A^'u^ 1 =l$~ l Rn +a5 _1 Ik 


uj , and A k_l u k_! = b k 1 


defines the equation at points on the coarse grid which have not been refined. Here R„ represents the global residua! after 
n sweeps on level k in the multigrid hierarchy. This indicates the level of convergence of the pressure, or velocity in 

Ak-l 

reversed flow regions, fields and is evaluated as R„ = b k - A^u^ . l£ 1 and Ik represent fine to coarse grid transfer 


Ak-l 

operators from level k to k-l. Ik is a simple injection operator that takes the solution from fine to coarse grids. In the 
current work, the covariant momentum balances in the (^,T|,0 directions are solved. There are two choices for the 

evaluation of on the coarse grids; (i) residuals of the covariant r\ and C, momentum balances can be directly 

restricted to the coarse grid, or (ii) the values of the residuals in the Cartesian balances can be computed through the inverse 
transformation. These values can be restricted to the coarse grid and then the residuals in the covariant balances can be 
recomputed on the coarse grid using coarse grid metrics. Numerical tests indicate that the second approach provides a 
much better representation of the fine grid problem on the coarse grid. This leads to significant improvement in the 
performance of the multigrid algorithm. Figure 3 shows a comparison of the convergence rates attained by the two 
approaches with the convergence history of a single grid computation (Curve 1). Clearly the second approach, restriction 
of the Cartesian residuals (Curve 3), is better than restriction of the covariant residuals (Curve 2). The example shown in 
Figure 3 was obtained for two dimensional flow past a trough geometry; similar behavior has been observed for all 
problems considered herein. 
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Since there can be sub domains within which the grid is refined in one or more coordinate directions with respect to their 
parents, it is necessary to implement all possible modes of multigrid transfer operations. For two-dimensional calculations, 
three different multigrid modes are implemented; namely, ^r| refinement, | refinement, and t| refinement. For three- 
dimensional SMGDD calculations seven different modes are implemented, i.e., T|, H,t|, £&, and refinement. For 

each mode appropriate interpolation schemes are implemented in order to provide initial guess for the fine grid based on 
coarse grid solutions. 



Figure 3. Comparison of convergence rates with different residual restriction. 

Instead of the standard V or W cycle (Brandt Ref. 14) a dynamic criterion is adopted in the present 
calculations to switch between the fine and coarse grids. When the ratio of the global residual between two successive 
sweeps is greater than a certain preset tolerance (typically 0.9 to 0.95), indicating that the convergence rate of the relaxation 
procedure is slowing down, the calculation is switched to the coarser grid level. If the ratio is greater than unity then the 
calculation continues on the fine grid since this could represent some flow feature (separation region, shock wave) is not 
completely settled in the fine grid computation. Hence it is necessary to allow a few more global sweeps in the fine grid to 
establish the solution (not completely converge) before switching to the coarse grid. Rubin et. al. (Ref. 3) have estimated 
the convergence rate of the global pressure relaxation procedure through a von Nuemann analysis of a linearized form of 
the RNS system of equations. The convergence rate X of the global relaxation procedure is determined by the spectral 
radius of the linearized system and is given by 

A.~ 1 -C,7t 2 (A^) 4 N § 2 /7i M 4 (6) 

where C, is a constant of 0(1), is the number of stations in the streamwise or marching direction, r| M is the extent of the 
normal boundary location and is grid size in the marching direction This analysis indicates that improved 
convergence is achieved when number of stations in the marching direction is reduced and when the normal extent of the 
compuiational boundary is small. Both these features are achieved in the present adaptive refinement calculations since the 
finer grids are limited to smaller regions of the flow field. Thus convergence rates on finer grids are comparable to coarse 
grids in most cases. The range of values for the residual ratio(0.9-0.95) has been established through numerical 
experimentation on a wide range of problems, although success of the method does not strongly depend on the choice ot 
this parameter. 

On the coarse grid it is not necessary to completely converge the solution. The tolerance for convergence on the 
coarse grid is prescribed to be an order of magnitude lower than the current norm of the residual on the immediately finer 
grid. The calculation also switches back to the coarser grid level if the magnitude of the correction to local boundary 
conditions on the fine grid is larger than the convergence criterion on the finest grid. 
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Figure 4 Comparison of eddy viscosity contours with and without local refinement. 

Applicability of the multigrid method to turbulent flows has been studied by many researchers, although most of 
the work has been limited to two dimensional flows. Braaten et. al. (Ref. 26) have applied the multigrid method to a 
pressure correction based solution procedure. Joshi et. al. (Ref. 27) have studied internal flows. Lein et. al. (Ref. 28) have 
also studied turbulent flows using the muitigrid method. Treatment of eddy viscosity for turbulent flows requires special 
attention. A significant difficulty with the multigrid method in turbulent flows, is the loss of resolution near solid walls on 
coarse grids especially when wall functions are not applied. This can result in erroneous values of eddy viscosity on 
coarser grids. Shyy et. al. (Ref. 24) have suggested a partial coarsening approach wherein grid points near the wall are not 
coarsened. This is a feasible approach but is difficult to implement in an adaptive refinement calculation as presented here. 
To the author’s knowledge the present study is the first that combines elements of adaptive refinement and the multigrid 
method to solve turbulent flow computations. Some new strategies to treat the eddy viscosity were tested and have been 
successfully applied to the turbulent flow cases presented herein. 

For all the turbulent computations presented herein, turbulence closure is provided through a Renormalization 
Group (RNG) based algebraic eddy viscosity model. The original RNG based algebraic model as derived by Yakhot and 
Orszag (Ref. 29) leads to a cubic equation for the eddy viscosity. This equation produces multiple positive roots in some 
cases, and the choice of the correct root is not clear. Lund (Ref. 30) has recast this cubic equation as a quartic that leads to 
four roots, two imaginary, one real positive and one real negative. Thus the choice of the correct root is clear. This form of 
the RNG algebraic model is applied herein. Details about the implementation of this model can be found in Kirtley et. al. 
(Ref. 31) 

In the present study, the eddy viscosity field is computed only on the finest grid and in those portions of the coarse 
grid where the grid is not refined. At coarse grid points that are refined, the eddy viscosity field is simply injected from the 
fine grid and held fixed during visits to the coarse grid. This strategy has been validated for flow in a two dimensional 
straight duct at Re=10\ Figure 4 shows the comparison of eddy viscosity contours obtained with and without local 
refinement. The agreement is excellent. The top figure corresponds to the adaptive calculation. The fine grid levels are 
limited to two narrow regions near the two walls. These are demarcated by the two horizontal lines. Note that the eddy 
viscosity contours are perfectly smooth across these boundaries. Before defining these thin regions near the wall boundaries 
the coarse grid solutions are obtained. On a coarse grid with no points in the laminar sub layer, the peak eddy viscosity 
values for this Reynolds number is computed to be about 200. The peak occurs roughly midway between the walls and the 
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centerline of the duct, a region of the flow field where no grid refinement is prescribed. On prescription of the thin 
refinement region and following the strategy explained above, it was found that the eddy viscosity values improves even in 
the region wherein no refinement is specified. On convergence, the eddy viscosity value and its location are captured 
correctly even though this location lies in the unrefined coarse grid outside the thin boundary layer region. This 
methodology has been applied successfully in all of the turbulent flow computations presented herein. In some 
computations, however, when grids normal to the surface are too coarse, the calculation failed to converge. In such 
situations, the grid normal to the wall is coarsened only once and then held fixed for all coarser grids and the grid 
refinement/coarsening is prescribed only for the streamwise direction. 


3. Governing Equations and Boundary conditions 

The grid adaptation methodology is demonstrated on the reduced Navier Stokes (RNS) set of equations. Detailed 
description of the derivation of these equations, associated discretization and applicability can be found in Ref. 2. One 
aspect of the RNS discretization is the use of consistent boundary conditions. The scheme is fully conservative because the 
discrete continuity equation is integrated directly and exactly at each grid point unlike pressure correction based methods 
that solve the Poisson equation for pressure correction to indirectly ensure mass conservation. This aspect provides 
robustness to the current adaptive gridding approach by ensuring mass conservation on local and global grids without the 
need for elaborate conservative interpolation schemes (Ref. 32, 33) to transfer boundary conditions from global coarse grids 
to local fine grids. Some elements of the boundary conditions are described here for completeness. 

Although, only internal flows are considered herein, boundary conditions will be described for both external and 
internal flows. There are situations where sub domains within the channel/duct may have a free boundary, and ‘external’ 
flow like boundary conditions are required in these local sub domains. The required inflow boundary conditions reflect 
whether the flow is compressible or incompressible. The free stream Mach number is prescribed for compressible flows 
and is set to zero for incompressible flow. 

At the inflow, U, V, W are prescribed for incompressible flow, (p is not required, but is obtained as part of the 
solution) and U, V, W and pgg = 0 for subsonic compressible flow. This allows for an adjustment in the inflow mass 


(through the density), that is consistent with the prescribed downstream pressure. In discrete form this equation becomes 


n 

P2 


Pi 


_n~l n 
P3 ~P2 


All A| 2 <6> 

where P 2 denotes the calculated pressure at first station downstream of the inlet and pj denotes the inlet pressure. The 
inflow pressure is given by p jn = (0 P| + (1 -Gd) p 2 . Here n denotes the global iteration counter. Note that the stations 
I and 2 are both at the current iteration level. A more detailed discussion of this condition can be found in (Rosenbaum et. 
a!.. Ref. 4). For supersonic flows the inlet pressure is treated as known since there is no upstream influence. For both 
compressible and incompressible flows the remaining boundary conditions are: 

U, V, W at all solid boundaries: 

U, W, p at all free (rp boundaries: ( V computed) 

U, V, p at all free (Q boundaries; ( W computed) 

At the outflow. (£=£; m ax) P 33 constant is prescribed; at an interface outflow, £=£*, < ^m ax , p is prescribed, either from the 
current grid level computation through inter domain transfers, or from multigrid transfers. This, depends on whether or not 
the boundary is shared with another sub domain in the current multigrid level. U, V, W are not required at any outflow 
boundary when U > 0; when U < 0 (this can occur if the sub domain boundary is located within a reversed flow zone), U, 
V, W are prescribed either from current grid calculation or from multigrid transfers. If the global domain ends within a 
reversed flow zone, then the negative convective fluxes are neglected. This is used to remove the need for a downstream 
velocity boundary condition necessitated by the upwinded streamwise convective terms. This has been validated for 
several flows with reverse flows at the outflow boundary. It is significant that the pressure at solid boundaries, V at free rj 
boundaries, and W at free C, boundaries are all calculated as part of the solution in a given grid level. For external flows, if 
a sub domain has its outflow at some (% < | max ) , then the boundary condition on pressure changes from Neumann to 


Dirichlet type. For internal flow, the outflow boundary condition for the pressure is always of Dirichlet type. In time- 
dependent, characteristic-based, Navier-Stokes computations that use locally embedded grids, boundary conditions are 
required for all variables, i.e., U. V, W and p on all boundaries. In addition, special care has to be taken to ensure that mass 
conservation is not violated locally or globally. With the pressure-based flux-vector splitting and the trapezoidal or 'box' 
formulation, this difficulty does not occur as the normal velocities V in T|, W in ^ or U in £*, are not prescribed at the (r|) 
upper or (Q cross flow, or (£,) outflow boundaries, but are computed as part of the solution. This eliminates the need for 
special interpolation formulae to ensure conservation of mass when the boundary conditions are prescribed from the coarse 
grid solution. Thus weak instabilities, that arise when such methods are applied in Navier-Stokes formulations without 
satisfying mass conservation, do not appear in the present method. Direct evaluation of the pressure at the inflow and wall 
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boundaries, with the trapezoidal or box discretization, also eliminates the need for extrapolated pressure boundary 
conditions at walls and outflow boundary conditions and ensures momentum conservation throughout. 

4. Results 

Several two and three dimensional flow problems, exhibiting incompressible to supersonic conditions, for both 
laminar and turbulent flow regimes have been considered. The choice of geometries is driven by availability of 
experimental data and/or calculations by other researchers. Some geometries have been chosen so as to bring out the 
efficacy of the segmented multigrid domain decomposition procedure and the ability to automatically define regions of high 
gradients that require improved resolution. As part of the segmented multigrid domain decomposition procedure presented 
herein, it is necessary to obtain solutions on a hierarchy of grids. The results on different grid levels can then be used as a 
means of demonstrating grid convergence. 

Unless specifically noted, all calculations, are initialized with uniform flow conditions on coarse grids. On the 
finer grids, interpolated solutions from coarser grids are used as the initial approximation to the solution. This is a 
significant feature, since in many full Navier-Stokes computations initial approximations must be obtained by Reynolds 
number continuation procedures (Ref. 34). This can add significantly to computational cost. 

4.1 Flow in an S-shaped duct 

Laminar and turbulent flow calculations are considered for an S-shaped duct as defined by Michelassi et. al. (Ref. 
35). Experimental data is available for both flow regimes. This geometry is typical of jet engine intakes that are required to 
redirect the flow. Secondary flow is a determining factor affecting the duct performance. Experimental data is available in 
Taylor et. al. (Ref. 36) for a Re=790, based on inflow hydraulic radius and bulk velocity. A partially developed inflow 
profile is assumed at the inflow. Figure 5 shows the three dimensional adapted grid for this geometry. There are four 
multigrid levels, of which two are adapted. The equivalent finest level corresponds to a 129x41x41 grid or 216849 mesh 
points. The adaptivity leads to three sub domains in the finest level with grids of 17x15x41 , 33x21 x41 and 17x19x41. This 
results in a reduction of nearly 75 %. The entire adaptive calculation (inclusive of all four levels) is completed in 7.45 hours 
on an IBM RS 6000-530 series machine or about one hour on a Cray Y-MP. Time under-relaxation was required for this 
geometry. A constant At = 1.0 was used on all grid levels. 

Figure 6 depicts the symmetric xy-plane grid, and the corresponding 'normal' velocity contours. Note that the 
contour lines are smooth across the sub domain boundaries. The local refinement pattern is defined by examination of the 
pressure contours that are shown on Figure 7. All the fine grid regions are automatically determined around regions where 
the pressure gradient is large. The smoothness of the solution across the domain boundaries is evident. Figure 8 depicts the 
cross-flow velocity vectors at various axial locations. The horizontal lines represent the boundaries of the local subdomains 
at these sections. The behavior of the cross flow changes significantly from section 2, to 3 and from 4 to 5. Note the 
presence of three counter rotating vortices in section 5. These are absent in the other three sections. Section 5 is located at 
the end of the second bend of duct. Figure 9 shows the comparison of the axial velocity component with the experimental 
data. The comparison is for all five sections at the z=0.5 or mid-span location. A grid convergence study is also presented 
on the figure. Good agreement is obtained with the data. The addition of fourth order damping was not required for this 
geometry. Also solutions were obtained with the RNS system without the comer region viscous terms in spite of the strong 
cross-flow separation that is present in this flow field. 

Turbulent flow calculations at Re=4().(XX) were also performed for this geometry. Figure 10 shows the velocity 
vectors along the symmetric xy-plane and the cross flow pattern at the end of the second bend (Section 5). The cross flow 
behavior is significantly different from the corresponding laminar flow, although recirculation regions are still present. At 
the inflow, a partly developed turbulent velocity profile is prescribed to match the experimental setup. Figure 1 1 shows a 
comparison of profiles of the axial velocity component at various sections along the mid span location with experimental 
data. Good agreement is again obtained. The turbulent flow computation was performed with two multigrid levels having 
17x21x21 points and 33x41x41 points, respectively. The eddy viscosity is updated globally after each sweep. This 
calculation was completed in about 7 hours on an IBM RS/6000 520 series machine or less than an hour on the CRAY Y- 
MP. Adaptivity was prescribed for a third level only in the streamwise direction. 

4.2 Flow in a 90° duct 

Ducts with rectangular/square cross-sections are very frequently found in engineering application. The laminar 
flow in a 90° bend duct, as defined by Humphrey et. al. (Ref. 37), is considered. The Re=790, based on inflow bulk 

velocity and hydraulic diameter (d) of the cross-section of the duct, and the Dean number De = Re 

where R c is the mean radius of curvature. 

A fully developed straight flow solution is prescribed as the inflow condition in order to match the experimental 
data of Humphrey et. al. (Ref.37). Figure 12 shows the three dimensional adapted grid for this geometry. Four multigrid 
levels are prescribed, two of which are adaptive. The finest grid is equivalent to a 65x33x33 mesh, although the adaptivity 
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actually reduces this number of grid points by about 60 %. The calculation is completed in 4.5 hours on an IBM RS/6000 
530 series machine or 35 minutes on a CRAY Y-MP. Note that the fine grid zones differ from those of the example shown 
in Figure 2. This is due to the difference in the inflow profiles. Figure 2 was obtained with a uniform inflow profile. As a 
result, the flow undergoes a significant gradient near the inflow; whereas, in the present case, the changes in the solution are 
negligible until the bend and maximum in the latter half of the bend. The adapted fine grid region is packed around this 
high gradient portion of the flow field. 

Figure 13 shows the contours of the streamwise velocity component (Vq) at four different cross-sections, viz., at 
the beginning of the bend (0=0°), two stations on the bend. (0=30°, 0=60°) and one at the end of the bend (0=90°). The x- 
axis is the normalized radial distance and the y-axis is the azimuthal direction. Note the change in the flow pattern. Figure 
14 shows the cross flow velocity vectors. The components of the velocity are the radial velocity (V r ) and azimuthal velocity 
w. Strong secondary recirculation develops in the second half of the bend. Such behavior has been observed 
experimentally by Humphrey et. al. (Ref. 37). The fluid tends to accumulate near the outside of the bend. Figure 15 
presents the comparison of streamwise velocity profiles at the four designated sections. The solution on two different grids 
is shown. Good agreement is obtained with the experimental data from Ref. 37. 

4.3 Three dimensional convergent divergent nozzle 

As a compressible flow example, the flow in a convergent divergent nozzle is studied for laminar and turbulent 
flow conditions. This geometry provides a good test case, as it spans the low subsonic to supersonic Mach range. The 
geometry has a square cross-section, as all four wails converge/diverge uniformly. For the present calculations, the ratio of 
the maximum to throat area is 3.81. At the outflow, a back pressure 0.04 times the inflow pressure is prescribed. This 
corresponds to the fully expanded case, as given by inviscid quasi one-dimensional analysis. 

For the laminar case, a Re=500 based on the throat hydraulic radius is considered. Figure 16 shows the adapted 
grid and velocity vectors on the symmetric xy-plane. Note that the refinement was concentrated around the throat region 
which represents the region of largest pressure gradient. Five multigrid levels were used. Two adaptive levels were 
specified. Although adaptivity was specified only in the streamwise direction, during the multigrid cycling the grid was 
coarsened in all three directions. The finest grid is equivalent to a 129x33x33 grid. The adaptivity reduces this number of 
grid points by nearly 55 %. The complete subsonic to supersonic, five multigrid, calculation is completed in about 14 hours 
on the IBM RS/6000 530, or less than two hours on the CRAY. Figure 17 shows the contours of the normal and streamwise 
velocity components. The solution is, once again, smooth across the zonal boundaries. The streamwise velocity increases 
nearly eight times from inflow to outflow. Figure 18 shows the contours of the normal velocity component at the throat and 
at the outflow. These flow patterns are quite different. This is also seen in Figure 19 for the cross flow velocity vectors. 
Note the presence of two counter-rotating vortices at the four corners of the throat section. Near the nozzle exit all the 
velocity vectors rotate outward as expected. Figure 20 shows a comparison of the streamwise and normal velocity 
components at the mid span (z=0) throat location for three different grids. No significant change is observed in the 
streamwise velocity component when the cross plane grid is refined, but axial refinement increases the peak velocity 
significantly. Note that the maximum velocity occurs away from the centerline. 

Calculations were also made for a turbulent Reynolds number Re=l() , and with several different values for the 
back pressure. Calculations are first completed for the fully expanded case, wherein Pback/Pin = 0-04. Figure 21 shows the 
pressure and Mach number contours on the symmetric xy-cross plane. Note that the flow smoothly accelerates from 
subsonic to supersonic Mach numbers. A magnification of the cross-flow near the four corners of the cross plane at the 
throat is also depicted. Two counter-rotating vortices are evident. This calculation was continued using three semi- 
coarsening multigrid levels (all full refinement), with a finest grid of 33x41x41 . Figure 22 presents more details of the flow. 
The top picture shows the velocity vectors in the symmetric xy-plane. Note that the flow accelerates significantly from 
inflow to outflow. The normal velocity contours are shown for three different cross-sections, viz., a section ahead of the 
throat, the section at the throat and a section near the outflow. Note that the flow behavior is quite different at the throat, 
wherein the recirculatory flow pattern as depicted in Figure 21 appears. This behavior is absent in cross-sections ahead and 
behind the throat. Although the contour line pattern for the sections in the convergent and divergent portion appear to be 
similar, the sign of v is reversed. The negative v contours are shown with dotted lines. The flow moves towards the 
centerline in the convergent portion and away from it in the divergent portion. 

4.4 Flow in an inlet geometry 

Computational fluid dynamic simulation of flow through engine inlets provides a testing methodology to analyze 
the performance of various engine configurations. The nature of the flow exiting from the engine will determine the 
performance of any turbomachinerv. 

A generic two dimensional inlet configuration has been selected. This consists of a 7.5° ramp on the lower 
surface, and a flat upper surface. The length of the passage is five times the inlet entrance height. The ramp terminates at a 
unit normalized distance from the inlet entrance. The calculation is for a free stream Mach number of 2.5 and a laminar 
Reynolds number of !0 3 ’ based on the height at inflow. The performance of the segmented multigrid domain 
decomposition procedure is assessed for this supersonic flow. Figure 23 depicts the adapted grid that corresponds to a 
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uniform fine grid of 161x81. Note that the fine grid is ‘packed’ in the regions where a shock is present. There are two 
regions of separated flow on the upper and lower walls that result from strong shock-boundary layer interaction. These 
regions are evident in the streamline pattern of Figure 24. The finest region in the adapted grid encompasses the two 
recirculation zones. Figure 25 compares the pressure contours obtained from the adapted calculation and full refinement 
calculation with five fully refined multigrid levels. Good agreement is obtained. The number of grid points for the local 
refinement calculation is about 8500. This compares with 13000 for the full refinement calculation. Due to the presence 
of the shock wave and associated reflections the reduction is less than in the earlier low speed examples. Since only 
‘rectangular’ regions have been defined as sub domain, the oblique nature of the shock wave leads to excessive refinement. 
This can be avoided by the use of an explicit solver (Ref. 13). this would allow each individual grid cell to be refined 
independently, and on the basis of a local truncation error criterion 

5. Summary 

A solution procedure that ensures accurate, grid converged flow solutions very efficiently has been described. It 
combines the advantages of the multigrid method and the flexibility of a domain decomposition approach to reduce 
computational costs through grid optimization. The overhead involved in computing truncation error and prescribing the 
finer levels of the multigrid hierarchy are minimal, leading to significant reduction in computational effort as compared to a 
non-adapted calculation with multigrid acceleration. Examples have been presented in incompressible and supersonic 
Mach number range. The pressure flux splitting based discretization scheme ensures that mass conservation is preserved in 
all subdomains. This is significant for the current procedure since boundary conditions have to be specified for local 
regions wherein a finer grid is prescribed and the need for special conservative interpolation schemes to calculate these 
boundary values is removed. The method also provided flexibility in choosing the adaptivity criterion. Higher order 
discretization can be applied in local regions of large gradients which are automatically identified in the present method. 
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Figure 5. Adapted grid for S-shaped duct; Four multigrid levels; Finest grid equivalent to 
129x41x41 grid points; Incompressible laminar flow at Re=790 



Figure 6. Adapted mesh and contours of v; Incompressible laminar flow in S-shaped duct 
Re=790 
























































Figure 15. Comparison of streamwise velocity profiles at different streamwise locations 
(midspan); Flow in a 90° bend duct; Re=790 
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Figure 16. Adapted mesh and velocity vectors on symmetric xy-plane; Compressible flow in 
a convergent-divergent nozzle; = 0.2, Re=500 





flow in a convergent-divergent nozzle; =0.2, Re=500 



Figure 18. Contours of v velocity component on cross-sections near the outflow and at the 
throat; Compressible flow in a convergent-divergent nozzle; =0.2, Re=500 





Three dimensional convergent-divergent geometry 
Crossections at the throat and near thd outflow 
Laminar flow, Re=500 

Solution on 129x33x33 grid (five multigrid levels) 
Cross flow velocity vectors 
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Cross-flow velocity vectors on cross-sections near the outflow and at the throat; 
Compressible flow in a convergent-divergent nozzle; = 0.2, Re=500 


Three dimensional convergent-divergent nozzle 
Laminar flow (Re=500) 

Grid convergence study 


Fine grid (129x33x33) 
Coarse grid 1 (33x33x33) 
Coarse grid 2 (33x 1 7x I 75 
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Figure 21. Pressure and Mach contours on symmetric xy plane, cross flow velocity vectors 
and streamline pattern near the comer for the cross-section at the throat; 
Compressible flow in a convergent-divergent nozzle; =0.2, Re=10 5 



cross-sections; Turbulent flow in a convergent-divergent nozzle; = 0.2 , 
Re=10 5 










Two dimensional inlet geometry 

Adapted Mesh (Five multigrid levels - Equivalent to 161x81) 








